This guided example is a supplement to the more general chapter offering a framework and workflow for power analyses. It follows the general steps outlined there. In this example, we will use the package brms for Bayesian estimation of parameter values. To keep things computationally tractable, we will focus on a simple comparison between two independent samples. The general steps are applicable to other Bayesian designs.
Pre-requisite: You will need to set up RStan on your computer for the Bayesian analyses in this example to work. See mc-stan.org and Rstan’s wiki.
I never managed to install “mac -os Rtoolchain” but I did manage to get #Rstan up and running thanks to the Coatless Professor’s workaround.
We will run a comparison of two independent samples, and seek to assess our capacity to detect 1 a range of possible effect sizes, considering Cohen’s \(d\) at 0, .15, .3.
#We will also increase the memory limit in R to allow for the large objects that are likely to be created:
memory.limit(100000)[1] Inf#This is a Windows specific function. Some possible ways to do this in Mac are suggested here: https://stackoverflow.com/questions/51295402/r-on-macos-error-vector-memory-exhausted-limit-reachedt.data.maker below
group.size (easy to adjust this to allow 2 sizes)effect.sizesim is just the counter of ‘which simulation this is’#| label: data-making-function
#| code-summary: "function to make simulated data"
t.data.maker <- function(sim, effect.size, group.size) {
# create a data frame with two groups, standard normal distribution
# group b differs from group a according to the specified effect size
data <- tibble(
a.group = rnorm(n = group.size, mean = 0, sd = 1), #create `group.size` rows of a standard normal,
b.group = rnorm(n = group.size, mean = effect.size, sd = 1), #... of a std normal shifted by `effect` for 'treatment group' (or after-treatment if this is within subject)
n = seq(1:group.size), #counter
'effect.size' = effect.size,
'sim' = sim
) %>%
pivot_longer(cols = a.group:b.group, #make this one row per 'observation in group' (rather than 1 per observation with 'both groups') ... so 2*group.size rows
names_to = 'group',
values_to = 'response')
return(data)
}We will run 250 analyses2 for each of the three effect sizes of interest. We use maximum group sizes of 1500.3 We put these into a format to pass to a map function, and then run the data-making function over the different effect sizes to generate the many data sets.
In the code below, we create this simulated data, creating a list of tibbles, one tibble for each simulated data set, considering the three effect sizes.
# set up furrr to run iterations more quickly ... this is (optional) parallel processing stuff
library(furrr)
plan(multisession)
options <- furrr_options(seed = 1234)
# make a tibble with our intended effect sizes (DR: for 250 simulations each)
effect.sizes <- tibble(
effect.size = c(rep(0, 250), rep(.15, 250), rep(.3, 250)),
sim = seq(1:750)
)
# run t.data.maker over all the requested data sets
# Map the list of effect sizes (and sim counter) to `t.data.maker`
#... to generate a list of tibbles (each of size 1500)
#... as specified, with 250 'treated' and 250 'controls' per effect size, where controls are standard normal and treated have the given effect sizes
sim.data <- future_map2_dfr(
.x = effect.sizes$sim,
.y = effect.sizes$effect.size,
.f = t.data.maker,
group.size = 1500,
.options = options
)
# split the simulated data sets (into a list of tibbles, one tibble per simulation) for analysis
sim.data <- sim.data %>% group_by(sim) %>% group_split()We set up a Bayesian regression model to run over the many simulated data sets.
It is important to ensure that the estimands we are interested in can be returned at this step, because this process can be time consuming. We don’t want to have to repeated it due to small mistakes. Hence, first carefully check that the code returns all the parameters you need, and then put it into a function that runs over the many data sets.
In any case, we need to run the model at least once outside of the function.4 This regression serves to compile the code and prevent the need to recompile it each time a new analysis is run, saving a lot of time.
Our simple formula for the t-test is a simple linear model:
reg.form <- response ~ 1 + groupIn the code below, we consider the nature of priors we can set for this model. We pass brms::get_prior our data, the family of response and link we want to use in the model, and the formula we are fitting. It suggests some possible priors or things to assign priors over?]
We aim for ‘weakly informative priors’, to constrain the model from entertaining highly unlikely parameter estimates. Setting reasonable priors on a model also helps speed up the MCMC sampling.
# load in requisite packages
library(brms)
library(tidybayes)
# check what priors can be defined in the model
get_prior(formula = reg.form,
data = sim.data[[501]],
family = gaussian()) prior class coef group resp dpar nlpar bound
(flat) b
(flat) b groupb.group
student_t(3, 0.2, 2.5) Intercept
student_t(3, 0, 2.5) sigma
source
default
(vectorized)
default
defaultNext we set our priors (coded).5
# set our own priors
(
reg.prior <- c(set_prior("normal(0 , 1)", class = "Intercept"),
set_prior("normal(0 , .8)", class = "b", coef = "groupb.group"),
set_prior("normal(1, 1)", class = "sigma"))
) prior class coef group resp dpar nlpar bound source
normal(0 , 1) Intercept user
normal(0 , .8) b groupb.group user
normal(1, 1) sigma userNote that we have specified
Now we want to check that the model runs properly, and that we can extract the estimands of interest. We run the regression with the brm() command, and we extract and manipulate the estimands of interest by pulling out the posterior distribution from the model with posterior_samples():
# grab a test data set from the many data sets
test.data <- sim.data[[501]]
# run the regression model
base.reg <- brm(formula = reg.form,
data = test.data,
family = gaussian(),
chains = 2,
cores = 2,
iter = 1200,
warmup = 200,
prior = reg.prior)
(base.reg) Family: gaussian
Links: mu = identity; sigma = identity
Formula: response ~ 1 + group
Data: test.data (Number of observations: 3000)
Draws: 2 chains, each with iter = 1200; warmup = 200; thin = 1;
total post-warmup draws = 2000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 0.01 0.03 -0.04 0.06 1.00 1790 1537
groupb.group 0.30 0.04 0.23 0.37 1.00 1804 1521
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 1.00 0.01 0.97 1.02 1.00 2239 1408
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).posterior <- as_tibble(posterior_samples(base.reg)) %>%
# I use the names of the posterior sample columns
# to calculate group means and Cohen's d
mutate(a.group = `b_Intercept`,
b.group = `b_Intercept` + `b_groupb.group`,
d = `b_groupb.group` / sigma)
(posterior)# A tibble: 2,000 × 7
b_Intercept b_groupb.group sigma lp__ a.group b.group d
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 -0.0316 0.318 1.00 -4260. -0.0316 0.287 0.317
2 0.0351 0.278 0.992 -4258. 0.0351 0.313 0.280
3 0.0170 0.317 0.990 -4258. 0.0170 0.334 0.321
4 0.0411 0.280 0.967 -4261. 0.0411 0.321 0.289
5 0.0241 0.311 0.965 -4262. 0.0241 0.335 0.322
6 0.0537 0.279 0.972 -4261. 0.0537 0.333 0.287
7 -0.00284 0.329 1.02 -4260. -0.00284 0.326 0.321
8 0.0180 0.317 0.999 -4258. 0.0180 0.335 0.318
9 0.0109 0.304 1.00 -4257. 0.0109 0.315 0.304
10 -0.0106 0.338 1.01 -4259. -0.0106 0.327 0.333
# … with 1,990 more rowsWe can see that our model run smoothly (decent effect sample sizes, Rhat of 1), and that we can return the parameters of interest. This approach can apply to essentially any Bayesian model - all one needs to do is know which part (or parts) of the posterior they are interested in and extract it. We can now move on to running the analysis over the many data sets. On my computer, this set of analyses took about 10 minutes to complete.
reg.maker <- function(simulated.data, base.reg = base.reg,
breaks = seq(from = 100, to = 1500, by = 100)) {
data <- simulated.data
# this function allows the sample size to be cut
# so that we can see power for many different sample sizes
cut.samples <- function(break.point, data) {
cut.data <- filter(data, n <= break.point)
cut.data <- mutate(cut.data,
sample.size = break.point)
return(cut.data)
}
data.cuts <- map_dfr(.x = breaks, .f = cut.samples, data = data)
data.cuts <- data.cuts %>% group_by(sample.size) %>% group_split()
# this function runs the regression model over the different sample sizes
run.reg <- function(data, base.reg) {
set.seed = as.numeric(data[1 , "sim"])
# note use of the 'update' function rather than brm
# to prevent need to recompile the model
output <- update(base.reg,
newdata = data,
chains = 1,
cores = 1,
iter = 1750,
warmup = 250)
output <- as_tibble(posterior_samples(output))
output <- mutate(output,
sample.size = nrow(data) / 2,
effect.size = pull(data[1, "effect.size"]),
sim = pull(data[1, "sim"]),
a.group = `b_Intercept`,
b.group = `b_Intercept` + `b_groupb.group`,
cohen.d = `b_groupb.group` / `sigma`)
return(output)
}
samples <- map_df(.x = data.cuts, .f = run.reg, base.reg = base.reg)
return(samples)
}
t1 <- Sys.time()
regression.samples <- future_map_dfr(.x = sim.data,
.f = reg.maker,
base.reg = base.reg,
.options = options)